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Abstract 

The event reconstruction in high energy neutrino telescopes is based on the expression 
of the time delay distribution of the detected photons propagating to some distance from 
their emission point through a uniform and scattering medium. 

Considering a realistic detector with a finite time resolution, we derive an expression 
for the probability density function of the time delay of detected photons. The asymptotic 
properties of this function as well as its corresponding cumulative probability distribution 
are calculated and discussed. 
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1 Introduction 



In a scattering medium, the Cherenkov photons emitted along the trajectory of a relativistic 
charged particle (or resulting from electromagnetic or hadronic showers of neutrino neutral and 
charged current interactions) and propagating to some distance are delayed relatively to the 
direct geometrical path. The expression for the time delay distribution of a photon reaching a 
detecting device is the essential ingredient in reconstructing events from high energy neutrino 
telescopes [1] and is used to calculate the likelihood of a given event hypothesis. 

In this letter, we consider a monochromatic and isotropic pointdike emitter-receiver problem: 
given a uniform medium characterized by scattering and absorbing optical properties and a 
gaussian time measurement uncertainty, we derive the arrival time delay distribution of the 
detected photons. The result is presented in a closed analytical form. 
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In the following section, we review the photon propagation function in a uniform medium [2] 
and show that it satisfies the expected limit of a zero time delay when there is neither scattering 
nor absorption, or when the emitter-receiver distance is zero. 

In section 3, we derive the probability density function (PDF) for a realistic setup, i.e. the 
detected photons are measured with a finite time resolution. We also discuss the properties of 
this function, useful in a numerical implementation. 

We conclude this letter with a discussion, presented in section 4. 



2 Photon propagation function p 

Scattering in a medium causes a photon arrival time delay t defined as a difference between the 
actual arrival time at the receiver, t arr i va i, and t geom , the arrival time from the path of a photon 
propagating in a medium without scattering: 

R . . 

t = t arrival tg eom t arrival ~ (1) 

where R is the distance between the emitter and the detection point, and c is the light speed in 
the medium. 

Our goal is to determine the PDF describing the distribution of time delay t of photons 
detected at the distance R from the emitter by the device with a finite time resolution. 

Instead of attempting the impregnable task of deriving the PDF from first principles, one 
might exploit the symmetry of the problem. This suggests the following equation for the photon 
propagation function p(R, i) [2]: 

p{R 1 + R 2 , t) = f dt'p(R 1 , t')p(R 2 , t - t') (2) 
Jo 

which describes the picture for p, analogous to the Huygens principle for waves in optics. Eq. (2) 
relies on the spherical symmetry of the medium and of the light source. p{R\ + R 2 , t) is the 
PDF corresponding to a photon propagation with a time delay t to a distance R\ + R 2 from the 
source (located at the origin). It can be expressed as the convolution of the PDF to propagate 
to the shell at distance R\ from the origin with the PDF to propagate from this shell up to the 
shell at distance Ri + R 2 from the origin. Once multiplied with an exponential damping factor, 
accounting for absorption and depending on the time delay, the normalized solution of Eq. (2) 
can be presented as 

P (e, p, t) = 1^5- I dt p& a o = 1 ( 3 ) 

where 

r R 1 c 

fs x- " S r + K < 4 ' 

and r(£) is the gamma function [3]. In (4), A a is the absorption length, A and r are phenomeno- 
logical parameters. When the photon propagation function p is derived from first principles, A 
and r are expressed in terms of the optical parameters of the medium and the differential cross 
section of the microscopic photon-medium interaction. As they stand in (3,4), these parameters 
must be extracted from the experimental data. 
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Clearly, for R — > 0, or for any distance without scattering, there should be no time delay, 
i.e. p must be concentrated near t — 0. These two regimes can be realized with the single 
requirement £ — > 0, provided that in (3,4) A is the effective scattering length in the medium 
(which is related to the mean free path and to the scattering profile in a way that when the 
effective scattering length increases, so does the mean free path). 

Examining the case £ — * shows that limg_ >0 p(£ ! P, t) does not exist - p is recognized as 
a generalized function with the kernel 7+(x) = £+ _1 /T(0 [4]. As is well known, generalized 
functions and their limits are not well defined for any values of their variables and parameters; 
the self-consistent approach is to consider (p, ip) = J pip instead of p, where ip(t) is a test 
function falling off at t — > ±oo fast enough, ensuring the finiteness of (p, ip) [4]. 

In terms of (p, ip), the condition t = for no scattering and/or zero distance can be readily 
verified: for any analytic function pit), we have 

lim (p, ip) = \imj o dtp(p, £, t) ip(t) = lim £ ^ <M£, P, *) *" = ¥>(0) (5) 

meaning that p has the limit (in the sense of generalized functions [4]) 

limp(£, p, t) = S + (t) (6) 

In (6), 5+(t) is a generalized function, filtering out the point t — +0 - origin on the t -axis, 
approached from the region t > (propagation time delay cannot be negative - p is defined on 
the right semi-axis t > 0). Eq. (6) supports our interpretation of A as the effective scattering 
length in a medium. 

The function p describes just the propagation in the medium, and can be interpreted as the 
PDF for the measured photon arrival time assuming an ideal detecting device. 



3 Photon arrival time PDF T a 

A PDF for realistic signals should account not only for the photon propagation and the proper- 
ties of the medium, but also for the features of the detector. Here, we consider the finite time 
resolution of the detector when measuring a signal. The function p described above can be inter- 
preted as a limit corresponding to an ideal measurement. The realistic PDF is thus constructed 
in terms of p and (p, the function simulating the detector time resolution distribution. 

We choose a gaussian distribution as a "detector function" p(t). The width a ^ of the 
gaussian distribution represents the time resolution and should be extracted from the character- 
istics of the detector, a is usually a composite from numerous sources affecting the final detector 
time resolution (which further motivates the use of a gaussian). 

Considering a detector with a finite time resolution instead of an ideal one is satisfactory 
from the physical point of view. We will see that this also removes the divergence at t — 0: the 
PDF accounting for a finite time resolution is an ordinary function which exists for all R, t. 

3.1 Definition of T a 

From the fact that a realistic setup will allow only the observation of the sum of the two in- 
dependent random variables, the photon time delay and the detector time resolution, it follows 
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Figure 1: p, emitter-receiver distance R = 45 m > R cr it = A (left), and p for R = 15 m < R cr i t = A, 
and 10 times p for R = 45 m (right). Divergence of p at small t and R < R cr it = A is shown. Throughout 
all the plots the values for the parameters are from [5]: A = 33.3 m, A a = 98 m and r = 557 ns. 



that the photon time residual PDF T a is the convolution of the photon propagation function p 
with the detector time jitter function g a : 

POO 

J^(£, p,t)= dt'p(£, p, t') g a (t - 1') (7) 
Jo 

where 

_ exp(-t 2 /2a 2 ) 
9 ° {t) ~ (8) 
incorporates the finite time resolution a^O. 

Since (7) is of the form (p, if), T a is not a generalized function, i.e. it exists for any R, t. 
It should be mentioned that when the time resolution is accounted for, the definition of the 
time residual t is slightly changed: it does not refer anymore to the photon arrival time but to the 
device trigger time - the arrival time corrected with a time jitter. Therefore t < is acceptable 
for J-'cr, in contrast to p, where the very definition of the time delay demands to consider only 
t > 0. 



x10" 3 




Figure 2: F a , emitter-receiver distance R = 15 m (left), and T a (solid line) and p (dashed line) for the 
distance R = 0.9 R cr it = 0.9 A = 30m (right). Time resolution a = 10 ns. Convolution is finite for any 
R > 0, — oo < t < oo, and the spurious singularity at R < A, t = 0, present in p, is absent in T a . 
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Now, the time residual t can also be negative, and the normalization condition 

/oo roo 
dtr a ( P ,z,t) = / dt' P ( P ,z,t') = 1 
-oo JO 



(9) 



is readily verified. 

From the physical point of view, negative time residuals do not contradict causality if T a 
decreases faster than g a , when t < 0. In other words, causality is preserved when g a envelops 
T a for negative t. In the next subsection, we show that this is the case. 

The integration in (7) can be carried out and we obtain: 



[no 'r ( irn 1 pt t 2 , ... 



(10) 



V^a 2 ~""V 4 2 
where D is the function of parabolic cylinder [3] and 

rj = pa — t/a 

For implementation purposes, it is convenient to express T a in terms of the confluent hyper- 
geometric function X F\ [3]: 



P, t) 



(ii) 



3.2 Properties and moments of T a 

The confluent hypergeometric function 1F1 is implemented in GSL, the GNU Scientific Library 
[6]. Here we cite a few useful properties of T a : 

• Large positive time residuals: t 3> a + per 2 . From the asymptotic for ±Fi [3], we obtain: 



= e^ 2 p(C,p,t) (1 



pa< 



5-i 



i + g2(1 -f-« (i- 



per 



2\ -2 



+ 



0" 



pa 



2 X -4N 



(12) 



• Behavior for negative t < 0. From the expansion of the function of parabolic cylinder D 
[3], we obtain for \t\ 3> cr: 



^(p, £, f) « (pa - ^)^exp(-t 2 /2 ( x 2 ) 



1 . ' po 



(13) 



Now, we are in the position to verify that causality is not violated: considering the ratio of T a 
to g a , we find that indeed T a decreases faster for negative time residual than the detector time 
jitter function: 

T a {t) ( _ t \~t 



9*{t) 



1 - 



per 



< 1 



(14) 
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• a < 1: using the relation lim CT ^o exp(— x 2 /2a 2 ) / V 2-rra 2 = 5(x), we recover p in the zero 
time resolution limit: 

lim^, p, t) =p(Z,p,t) (15) 

a— >0 

• T a is finite for t — and any R; e.g. for R = (£ = 0) 

^ CT (0, p, 0) = — L=; lim ^(e, p, = tf(f) (16) 



• Long distance: for £ ^> 1 we have: 

• Short distance: in the limit ^ < 1, the receiver time response function survives: 
exp(-t 2 /2(T 2 ) " 



(17) 



•?>(£, P, 



1 + |f 2 v ^ln(p ( r) + lE + 2 In 2 - 7rer/i(^)j + 0(£ 2 ) 



(18) 



where erfi(z) is an error function of the imaginary argument, erfi(z) = ^= J Z dxe x2 , and 
7s w 0.57721 is the Euler's constant [3]. 

These properties are useful to implement algorithms and to avoid numerical overflows. E.g., 
for large t, the calculation speed, as implemented in GSL [6], of \Fi(a, b; x) slows down with 
increasing x and numerical overflow may occur. This can be circumvented by switching to the 
asymptotic behavior of T a for large t. The case of large £ can be treated using (17) to avoid 
numerical overflow caused by T(£), appearing in (11) and in the expression for the cumulative 
of T a (see below). 

Also interesting are the moments of T a : 

/oo 
dtt N ^(p,Z,t) (19) 
-oo 

which can be calculated in a closed analytical form; we explicitly present here the first two - the 
time residual average and the corresponding dispersion: 

(*> = -, ((Atf)^(t 2 ) - (t) 2 = a 2 + ^ (20) 
P P 

As expected, only the detector time resolution survives in the variance at £ = 0. 
3.3 Cumulative of T a 

T a is a PDF for a single detected photon. Another object of interest is the time delay probability 
density of the first photon in a sequence of N detected photons. It is given by 

Fi N \p, tt)=N- Mp, e, t) ■ (1 - C^P, £, t)j (21) 
where Cjf is the cumulative probability: 

Cr(p,S,t)= f dbF {p,Z,t) (22) 

J — oo 
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Figure 3: Cumulative PDF Cjr versus t (ns) for R = 50 m. Time resolution a = 10 ns. 



This function appears in the expression for the probability u(p, £, a, t; e) that the time 
residual lies between t and t + e: 



o>(p, £, <r,t;e) = J t 1 dt'F a (p,^t') = Cr(p, £,t + e)- C^P, £, t) (23) 

Note that since JF CT is a smooth function defined everywhere, for e t the probability above can 
be approximated by T a : 

u(p,^a,t;e)^e^ a (p,^t + e/2) (24) 

This prescription cannot be used in the case a = 0. Instead, the analytical solution for the 
cumulative distribution C p , given in terms of incomplete gamma function [3] T(£, pt), must be 
used: 

c p = f dt' P (p, t, f) = i-r& pt)/m 

Jo 

The analytical expression for Cjr with an arbitrary £ cannot be obtained. From (11), it follows 
that the cumulative distribution can be expressed as an integral involving the error function: 

CA P ,C,t) + ^ f d yy^e^erf(^) (25) 



erf I -J=) - exp(pV/2 - pt) I 1 - erf[ P] 



When £ = 1 the result is 

CAM ' t] ^ = \ + \[-'\^) -'VV2 

The case of integer £ can be handled using 

roc f) n i T (in 

jf dyy n e-^erf(cy + b) = (-1)"^- er/(6) + exp((p 2 +46c)/4c 2 ) f 1 - er/(& 



(26) 



(27) 



An analytic continuation for a fractional order derivative is possible [4], but not practical 
as higher hypergeometric functions emerge. In the Appendix, we give a simple semi-numerical 
approximation valid for any £ > 0. 

In [5], further examples involving the cumulative distribution are presented. 
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4 Discussion 



In this letter, we suggested and calculated the arrival time probability density for a photon at 
some distance from an isotropic point-like source in a uniform scattering medium for a realistic 
detector. JF CT , accounting for a finite detector time resolution a ^ 0, is free from the spurious 
singularities which appear in the generalized function p (the unrealistic zero time resolution limit 
of To) and is defined for all t and R > 0. The expression for T a has been obtained in a closed 
analytic form (10, 11) and contains parameters characterizing the medium and the detector. We 
further discussed properties of useful for a numerical implementation. Finally, we calculated 
semi-analytically the corresponding cumulative probability entering multi-photon PDF's in order 
to provide a usable numerical implementation scheme. 
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Appendix: Numerical integration scheme for Cj= 

We sketch the procedure for the calculation of the cumulative Cjc. 

To avoid numerical integration from to oo in (25), we use the asymptotic properties of the 
error function [3], and split the integral on three subintervals: (0, oo) = (0, y m \ n ) © (y m in, l/max) © 
(l/max, oo). The parameter a in y m i n = max (0, t — a V2<r 2 ) and y max — max (0, t + a V2a^) is 
chosen from the condition ^7 y \ > en, leading to erf((t — y)/ V2a 2 ) m ±1. E.g., choosing a = 2.5 
results in an error less than 0.041%. 

The integration on the first and third subintervals, where erf — ±1, can be carried out and 
we are left with an integration over the finite interval (y m m, 2/max): 

CAP , e . t) - l - r«.^(^...) + J; Jm £■ d „ s . 

where T(a, x) = f£° dz z a ~ 1 e~ z , the incomplete gamma function [3], is implemented in GSL [6]. 

Integral with finite limits J is finite and can in principle be calculated numerically, but when 
£ < 1, numerical integration algorithms may become unstable when y m i n — * 0. This can be 
circumvented by handling J separately for the two cases: 

1. s > l. 

J is calculated numerically. Since y m i n is positive and £ > 1, there is no problem with the 
convergence; in this standard integration algorithm [7] can be applied. 

2. £ <1- 

When y min <^ 1, y^ is an integrable singularity in the region y ~ y m m- In order to isolate 
explicitly the terms ~ Z/min) a small finite e, y m \ n < e < y max , is introduced and the 
integration is performed analytically in the range y m i n < y < e. Expanding exp(— py) and 
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erf((t — y)/\/2a 2 ) in y, we obtain: 

■A 



J 



2r(0 



+ 



2r(0 



2/min 
!/max 



V2V 2 



y^ 1 e ra er/ 



+ 



+ L dyy '" (1 - m) { erf {^)- y ' 



-t 2 /2<7 2 



2r(0 

+ erf 



t yi 



V2 



7T(7 Z 



3 -t 2 /2<r 2 



+ 



+ per/ 1 



pHL + 0(^+ 2 ) 



and the remaining integral 

dyy^e-py erf((t - y)/V2o*) 
can be calculated numerically using standard integration algorithms. 
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